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A method is introduced for the construction of meshless discretization schemes which pre- 
serve Lie symmetries of the differential equations that these schemes approximate. The 
, method exploits the fact that equivariant moving frames provide a canonical way of associ- 

■ ating invariant functions to non-invariant functions. An invariant meshless approximation 
of a nonlinear diffusion equation is constructed. Comparative numerical tests with a non- 
invariant meshless scheme are presented. These tests yield that invariant meshless schemes 

\ can lead to substantially improved numerical solutions compared to numerical solutions 

■ generated by non-invariant meshless schemes. 

^ ; 1 Introduction 

Invariant discretization schemes have received increasing attention over the past 20 years, see 
e.g. [3,4,6,7,13,17]. Such schemes are attractive in that they preserve an important property 
of a system of differential equations, namely its maximal Lie invariance group G or at least a 
certain physically interesting subgroup of G. 

One of the most distinct properties of invariant discretization schemes for evolution equations 
is that they in general require the usage of a moving discretization mesh. The necessity of 
, using non-orthogonal and/or non-stationary meshes considerably complicates the construction 

and analysis of invariant numerical schemes, especially in the multi-dimensional case. Special 
techniques can be used to overcome this problem, such as the symmetry-preserving discretization 
in computational coordinates [4,11] or invariant interpolation schemes [3]. A problem with the 
\ former technique is that in the higher-dimensional case quite some computational overhead might 

be required to construct a proper mapping from the computational domain to the physical space 
of the system of differential equations. 

On the other hand, moving meshes or grids that are adapted to complicated domain geome- 
tries which hamper the straightforward use of finite differences or related discretization strategies 
are not new in the numerical analysis of differential equations. In fact, constructing, storing and 
modifying discretization grids is costly and hence is attributed as one of the main drawbacks 
of the otherwise popular finite element method [5, 15, 16]. Indeed, often a large fraction of the 
computational time spent for the numerical integration of a system of differential equations is 
consumed by the construction of the discretization mesh itself. 

This is why a new class of discretization schemes is steadily growing in importance over 
the past years, namely so-called meshless schemes. The main observations on which meshless 
methods rest is that no information on the connectivity of the single nodes at which the numerical 
solution is sought is required in order to discretize a system of differential equations [5, 15, 16]. 
All the information needed to approximate the various derivatives is already included in the 
nodes themselves. This makes meshless methods attractive to (at least partially) avoid the 
computational overhead required by the construction of discretization meshes. 

It is then obvious to ask whether meshless methods could be employed in the construction of 
invariant numerical discretization schemes as well. This would allow one to (partially) bypass the 
complication one faces in attempting to find discrete approximations of a system of differential 
equations that preserve the symmetries of that system. 



1 



It is the purpose of this paper to describe an algorithm for the construction of invariant 
meshless discretization schemes. The key idea on which our construction relies is a property 
of equivariant moving frames to send a given function to an invariant function [9, 10]. This 
property was successfully exploited to construct invariant finite difference schemes for partial 
differential equations starting from non-invariant schemes [3,12,21]. In the present paper we 
extend this method to meshless discretization schemes. 

As what concerns the organization of this paper, in Section 2 we describe the general proce- 
dure for finding invariant meshless discretization schemes. This method is applied in Section 3 
to construct an invariant meshless scheme for a nonlinear diffusion equation. Numerical tests 
comparing the invariant with the non-invariant meshless scheme are carried out in Section 4. 
The conclusion and some thoughts for further research directions are presented in Section 5. 



2 Invariant meshless discretization schemes 



There is not a unique way of constructing meshless approximations to differential equations. In 
fact, a number of different strategies to discretize a differential equation without or only partial 
usage of a discretization lattice are used, such as meshless (generalized) finite differences, smooth 
particle hydrodynamics, the element free Galerkin method or the meshless local Petrov-Galerkin 
method. For a review of these and further techniques, see e.g. [2,5, 14-16] and references therein. 

In the present paper we exclusively work with meshless discretizations based on meshless finite 
differences. We stress, though, that similar techniques as introduced below could be applied to 
other meshless methods that discretize a system of differential equations in the strong form. 

Meshless finite differences basically rest on the expansion of a function u:W—, 
multi-dimensional Taylor series around the node x°, 
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Truncating this series to the mth order derivatives one is left with an expansion that includes 
s = (p + m)\/(p\m\) coefficients. Thus, in theory s nodes x J are needed to solve for the s 
derivatives u a evaluated at the node x° from the linear system 
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j = 1, . . . , s. The s nodal points (x^ ,u 3 ) are usually chosen to be neighboring points lying inside 
the p-sphere of radius r centered around the node x . In particular, the nodes are not required 
to lie on a predefined, topologically connected mesh which makes the method truly meshless. 

The practical problem that can arise in this construction is that for certain distributions of 
the s nodes (e.g. all points lying on a line), the system (1) cannot be solved for the required 
derivatives u a \ x o. A possible ad hoc remedy for this problem is to include more than s points in 



the system (1), i.e. to over-determine it [5 
squares solution of (1), which reads 

( u ( m) ) d \ x o = (S T WS)~ 1 S T Wb, 



14]. The derivatives u a \ x o then follow from the least 
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The vector (u^) d \ x o contains the s derivatives of the truncated Taylor series evaluated at the 
node x° where Uj is the highest derivative occurring. The matrix S is build from the associated 
coefficients of these derivatives where Ax„ = (x- 7 — x°) a . The vector b includes the k > s 
functions values v? of the nodes inside the radius r of the p-sphere centered at x°. A (diagonal) 
weight matrix W is included in the least squares solution for the geometrical reason to give 
greater weight to the points x^ closer to x . More details on this construction can be found 
in [5, 14]. The extension to vector- valued functions u = (u±, . . . , u q ) is straightforward. 

The meshless approximated derivatives can be used to discretize a system of differential 
equations C: Ai(x,u^ m ^) = 0, I = 1,...,L, where u^ m ' includes all the derivatives of u with 
respect to x up to order m as well as u itself. This leads to a meshless numerical scheme 
V: Di(xi , (u^) d ) = 0, where {u^ m " > ) d denote the discretizations of derivatives in using 
Eq. (2). If the system C includes derivatives of u with respect to the time t then D in addition 
to the meshless spatial derivatives requires a discretizations of the time-dependent derivatives. 

We now briefly describe the method of invariantization using moving frames. An extended 
discussion can be found in several excellent papers on that subject, including [9, 10, 17, 18, 21]. 

Definition 1. A (right) moving frame p is a smooth map p: M — >• G from a manifold M to a 
finite dimensional Lie group G acting on M with the property that 

p(g ■ z) = p(z)g~ 1 , (3) 

for all z 6 M and g 6 G. 

The theorem on moving frame requires a group action to be free and regular in order to 
guarantee the existence of the moving frame p. It should be noted that in Definition 1 the group 
G is restricted to be finite dimensional. This restriction is in fact only apparent as the theory of 
moving frames is already formulated for infinite dimensional Lie (pseudo)groups, see e.g. [19]. 

If the group action of G is not free, it can be made free by constructing the moving frame on 
a jet space J m = J m (M, p) of appropriate order m. An alternative way of making a group action 
free is to extend it to the product action on the Cartesian product of copies of the original 
space M, denoted by M° = {(z 1 , . . . , z k )\ x l ^ x^ for alH ^ j} where z J = (x^,vP) are the 
nodal points. The joint product action is simply the component-wise action, g ■ (z l , . . . ,z k ) = 
(g ■ z 1 , . . . ,g ■ z h ), see also [3, 17, 21]. 

Moving frames are determined from a method referred to as normalization. For a r di- 
mensional group action, in this procedure one sets up a system of r equations involving the 
coordinate functions z, where z £ {xi,u^} in the case G acts on the mth order jet space J m 
or z € {x\,ui}, in case the action of G on the joint product space M° is considered. The first 
possibility leads to a moving frame 

p {m) Qn 

the jet space J m , p^ : J m — >• G while the second 
possibility leads to a product frame p° on the space M°, p° : M° — > G. 

The system of r normalization conditions is chosen in such a manner as to determine a sub- 
manifold of J m (or M°) which intersects the group orbits only once and transversally. One usu- 
ally sets r of the coordinate functions z to appropriately chosen constants, i.e. the normalization 
equations are z\ = c\, . . . , z r = c r , although in the discrete case it is beneficial to set combina- 
tions of the coordinate functions to constants, see the example in Section 3. Then one replaces 
these equations with their respective transformed forms, i.e. Z\ = g-z\ = c±, . . . , Z r = g-z r = c r , 
and solves this algebraic system for the group parameters g = (ei, . . . ,e r ) in terms of z. The 
result of this construction is the right moving frame p^ (or p°). 

For the present purpose the most important property of moving frames is that they define a 
canonical map from a given (non-invariant) function to a G-invariant function. 

Definition 2. The invariantization of a real-valued function / : M — > R using the (right) moving 
frame p is the function t(/), which is defined as i(f)(z) = f(p(z) ■ z). 
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Invariantization is the key for constructing invariant discretization schemes. To accomplish 
this, the product frame p° is computed using the symmetry group G of the system of differential 
equations C and extending its action to the joint space M° . The product frame should be 
compatible with the moving frame p^ m \ which is determined for the action of G prolonged to 



the jet space J m . Compatibility means that p° 



P 



(m) 



in the continuous limit {x J } 



x. 



This 



is achieved by computing the product frame p° using the discretized form of the normalization 
conditions that are used to construct the moving frame p^ m \ see [3,17,21] for more details. 

With the moving frame p° at hand one can invariantize any standard numerical scheme 
T> which approximates a system of differential equations C. The invariant scheme associated 
with T> is i(T>): Di(i(x^), i((u^ m ') d )) = 0, I = 1,...,L, which in the continuous limit {x^} — >• 
x yields the system of differential equations C expressed in terms of differential invariants, 
l(C): Ai(l{x),l(u^)) = 0. See [3,12,17,21] for further details. 

The extension of the invariantization procedure to meshless discretization schemes is now 
straightforward. Once the moving frame p° is determined on the space of nodes (x^u 3 ) it can 
be applied to the least squares solution (2) by invariantizing the vectors (u^ m ^) d \ 
the matrix S in the following way, 
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The invariantized derivatives in L((u^ m ^) d \ x o) and the invariantized nodal points (t(a^'), t(u 5 ')) 
are sufficient to invariantize any meshless approximation T> of the system C This is done by 
replacing the occurring nodal points and meshless discrete derivatives in T> by their respective 
invariantizations, just in the same way as this is done in the case of finite difference schemes. 
An example for this construction is presented in the subsequent section. 



3 Invariant meshless scheme for a nonlinear diffusion equation 

In this section an invariant meshless scheme is constructed for the nonlinear diffusion equation 
u t = (u-* /3 u x ) x . (4) 
The meshless Euler forward scheme for Eq. (4) is given by 

^ = -^- 7/3 (^) 2 +-- 4/3 -i, (5a) 

where u and u stand for the values of u = u° at the subsequent and previous time step of the 
integration, respectively, and At = t — t is the (constant) time step. The meshless derivatives 
u x and u xx are evaluated at t and x = x°. Likewise, the meshless leapfrog scheme for Eq. (4) 
reads 

^ = -\u-y\uif + u-^{ut). (5b) 

In the numerical results reported in Section 4 we use the leapfrog scheme (5b) and every 20 
steps the Euler scheme (5a) to suppress the computational mode of the leapfrog integration. 

In both these schemes the derivatives u x and the meshless approximations of the 

derivatives u x and u xx . These derivatives are found from the least squares solution (2), where 



4 



in the present one-dimensional case 
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As a weight matrix we use W = diag(exp(— /i(Ax J ) 2 /r 2 )), j = 1, . . . , k and \i = const. 

The reason for also including the third derivative in the vector (u^) d \ x is that it increases 



the order of approximation of the derivatives u x and u%. x . If the matrix S is square, it can be 
directly inferred from the truncated Taylor series that (at least on a uniform grid) including u xxx 
leads to first derivatives with third order accuracy and to second derivatives with second order 
accuracy. Moreover, in [5] it was shown (again for a uniform grid) that the least square solution 
(S non-square) invoked for finding {u lym ^) d \ x does not degrade the accuracy of the approximation. 

We now discuss the invariant meshless approximation of Eq. (4). The diffusion equation (4) 
admits a five dimensional maximal Lie invariance algebra q which is generated by 



d t , d x , 2td t + xd x , 2xd x -3ud u , 



x 2 d r 



3xud u , 



see e.g. [8] where an invariant finite difference discretization for (4) was constructed. The single 
vector fields in q can be exponentiated to one-parameter Lie groups, which can be composed 
to yield transformations from the five-parameter maximal Lie invariance group G of (4). The 
transformations of G acting on M = {(t,x,u)} are of the form 



T = e 2£3 (t + ei), X = e £3+2£4 + 

\l-e 5 x J 
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-3£4 I 
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£x,x) 3 U. 
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The action of G becomes free on the first jet space J 1 = J 1 (M, 2). The moving frame on J 1 
is constructed from the normalization conditions t = 0, x = 0, u = 1, ut = 1 and u x = 0. This 
allows one to solve for the group parameters s±, . . . ,£5. The resulting moving frame is 



x 2 u x + 3xu 




(7) 



xu x + 3u 



We now turn to the construction of the compatible product frame p°. The joint product 
action of G on M° follows from evaluating (6) at the single nodes (t 3 , x 3 , v? ) , i.e. 



T* = e e3 (t j + £1), X* 



3 £3+2£4 



1 - e 5 x 3 



+ e 2 
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£ 5 X j ) 3 U j , 



(8) 



and similar on the subsequent and previous time layers. It is readily verified that the joint 
product action leaves invariant the condition for the nodes to remain fixed during the integration, 
which is x 3 — x 3 =0. For this reason, fixed nodes do not break the invariance of Eq. (4) and 
we can assume that x 3 = x 3 , which we do for the sake of simplicity. The same is true for the 
equation t 3+l — t 3 = 0, which defines the flatness of the time layers. Hence, t 3+l = t 3 = t. 

Before constructing p° it is worthwhile pointing out that the scheme (5) is already invariant 
under the action of the one-parameter groups associated with ei, . . . , £4. This can be seen as both 
the variables t and x 3 only arise in the differences At and Ax 3 , which are obviously invariant 
under translations in t and x. At the same time, the scaling properties of the scheme (5) 
are the same as of (4), which implies the scale invariance of (5). Thus, the only symmetry 
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transformation which is violated by (5) is associated with the group parameter e 5 , i.e. inversions 
in x. This is why we construct the invariantization map only for this group parameter. 

In the moving frame p^ l \ the component £5 follows from the normalization condition u x = 0. 
This normalization condition is replaced by u d = to guarantee the compatibility of the moving 
frame p° with p^ in the continuous limit. The problem when using the least square solution (2) 
to obtain u d is that it will be very hard (or even impossible) to find the moving frame component 
£5 as the normalization procedure boils down to solving a high order polynomial equation for £5. 
This is why we use a less accurate approximation of u d to compute £5. In particular, using the 
Taylor series expansions u r = u + u d {x r — x) and u 



at x, where (x T ,u r ) and 



[x l , u l ) 



u — u d (x — x ) is sufficient to determine 
are the nodes lying immediately to the left and to the right 



of x. This leads to the usual centered difference approximation 



x' 



(9) 



Note that on a non-uniform grid it is not guaranteed that this approximation is second order 
accurate. However, it is numerically verified in the following section that the invariant meshless 
scheme which follows from the moving frame that employs the approximation (9) for u d is more 
accurate than the associated non-invariant scheme (5). 

Computing the p°-component £5 from the normalization (9) leads to the cubic equation 



u r (x r ) 3 — u l (x 1 ) 3 



u r (x r ) 2 — u l (x 1 )' 



d + 3- 
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(10) 



This equation is solved using the explicit formula for the roots which gives three distinct solu- 
tions. In the numerical tests we carried out only one of the roots was real, which was then used 
as the moving frame parameter £5. 

This apparent non-uniqueness of the component £5 in p° already arises when computing £5 
in the moving frame p^ from the normalization u x = 0, which in addition to the solution 
£5 = u x /{xu x + 3u) also has the solution £5 = 1/x. This solution for £5 cannot be used as then 
the normalizations x = 0, u = 1 and Ut = 1 could not be solved for £2, £3 and £4, respectively. 



Indeed, it is verified that in the continuous limit x r — x — > 0, x 



0, Eq. (10) has the root 



1/x of multiplicity 2 and u x /(xu x + 3u) as a simple root. This shows that the moving frames 
p° and p^ are indeed compatible. 



The invariant meshless counterpart to the scheme (5) is given by 
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The invariantized discrete derivatives follow from solving Eq. (2) using the invariantized expres- 
sions for (u^) d \ x , S and b, which are, respectively, 
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Likewise, the weight matrix W is invariantized to give t(W) = 
these formulas, £5 is the real solution of Eq. (10) and we have 



V t(n 3 ) ) 
diag(exp(-/it(Ax : ') 2 /r 2 )). In all 



i(x j ) 



1 - e 5 xi 



i{ui) = (1 



£ 5 x j ) 3 u j , 



t(u) = (1 — e^x) 6 u, l(u) = (1 — e^,x) 3 u. 
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4 Numerical tests 



We compare the invariant scheme (11) against the non-invariant meshless scheme (5) by carrying 
out numerical tests with the following three exact solutions of Eq. (4), 



Ul 
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(2cix — 3cft + C2) 
(x + c^ 2 



-3/4 



t + C 2 

u 3 = (cix + c 2 ) 



+ c 3 (t + c 2 ) 



-3/4 



(12) 



where ci, 02,03 are arbitrary constants. The third solution is a stationary solution. For these 
and further solutions of Eq. (4), see [1,20]. 

In all the numerical experiments reported we solve Eq. (4) with the invariant meshless 
scheme (11) and the non-invariant meshless scheme (5). We carry out the numerical integration 
on the interval L = [1,2] and choose the constants ci, 02 and C3 so that the respective exact 
solution u e is not singular within L. On this interval, we first create an equally-spaced grid with 
N = 40 grid points. Each of the grid point is then perturbed by adding a Gaussian distributed 
random number with zero mean and standard deviation 0.1 • Ax, where Ax is the spacing of 
the initial uniform grid. Dirichlet boundary conditions are used with the values of vP at the 
boundaries given by the corresponding value of the exact solution ui , U2 or u% . 

Ten independent integrations using ten different realizations of the above described grid 
generation procedure are carried out for 1000 time steps of the size At = 0.001. The time step 
is rather small so as to avoid numerical instability in the course of the integration. Larger time 
steps could be used if implicit schemes would be employed instead of explicit schemes in (11) 
and (5). The resulting root mean square errors 



rmse 



1 n 

- £K 



J\2 



are computed after each integration, where v? n and ul are the numerical and exact solutions of 
Eq. (4) at t = 1 at the nodal point x J , respectively. The averaged root mean square errors for 
the three solutions (12) for the non-invariant and the invariant scheme are reported in Table 1 
and denoted by rmse n j s and rmse^, respectively. 



Table 1. Root mean square errors for the exact solutions (12) of the nonlinear diffusion equation (4). 
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In the first run (solution ui), ci = 02 = 0.1, in the second run (solution U2), c\ = C3 = and 
02 = 10 and c\ = C2 = 0.1 in the third run (solution u$). It is worthwhile pointing out that 
the weight matrices W were chosen differently for the invariant and the non-invariant scheme. 
The reason for this is that the invariantization of Ax J enters the weight matrix of the invariant 
scheme. To facilitate the comparison of the results, we have tuned the parameter fi in each of 
the runs of the invariant scheme so that the entries in the weight matrices of both the invariant 
and the non-invariant scheme are of the same order of magnitude. 

Table 1 shows that the invariant scheme is able to better approximate the exact solution at 
t = 1 in all three test cases. On top of that, for the case of the stationary solution U3 we found 
that the invariant scheme approximates the exact solution up to machine precision. 
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As a further sensitivity test we run several integrations of the invariant scheme (11) and the 
non-invariant scheme (5), respectively, and vary the parameter \i in the weight matrices W. 
This parameter controls the influence of distant grid points in the meshless approximation of 
the discrete derivatives of u at the center node. We use u\ as the exact solution in these runs. 
In Fig. 1 we depict the result of this sensitivity study. As was discussed before, it is necessary to 
use different \i in the invariant and non-invariant integrations due to the different magnitudes 
of Ax J and t(Ai J ). It can be seen from Fig. 1 that the invariant numerical scheme (11) gives 
(substantially) better integration results over virtually all values of the parameter fi. 



His 

3 7 11 15 19 23 27 




2 4 6 8 10 12 14 

H„i, 



Figure 1. Sensitivity study of the invariant scheme (11) and the non-invariant scheme (5) with respect to the 
parameter in the weight matrix W . Upper x-axis: /j, used in the invariant scheme. Lower x-axis: fi used in the 
non-invariant scheme. 

The final sensitivity test we carry out is with respect to the parameter r, i.e. the radius within 
which grid points are used to compute the meshless approximations of the partial derivatives of u 
at the node x. The result of this study using u\ as the exact solution is depicted in Fig. 2. Again, 
varying r the invariant scheme is (substantially) better than the non- invariant. Moreover, for 
the lowest value of r chosen, r = 0.16, the non-invariant scheme (5) did not converge, whereas 
the invariant scheme (11) produced the best root mean square error of all the integrations. This 
means that fewer nodal points are needed to compute a stable approximation of the derivatives 
for the invariant scheme in the present example. 

5 Conclusion 

In this paper we developed a technique for the construction of invariant meshless discretization 
schemes. The key of this method is the application of the moving frame invariantization map 
to meshless discrete derivatives which ultimately boils down to the invariantization of a system 
of truncated Taylor series expansions. 

Despite practically demonstrated solely for a one-dimensional evolution equation, the method 
of invariant meshless discretization as introduced in this paper is particularly suitable for multi- 
dimensional problems. Differential equations with more than one space dimensions posed a se- 
vere problem for the invariant scheme construction machinery available so far due to the necessity 
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Figure 2. Sensitivity study of the invariant scheme (11) and the non-invariant scheme (5) with respect to the 
parameter r. 

of using non-orthogonal, possibly moving discretization meshes. Consequently, most of what is 
known by today about invariant discretization schemes has been learned from the consideration 
of ordinary differential equations or single (l+l)-dimensional evolution equations [6,8,12,13,21]. 
It was only recently that methods for the construction of invariant numerical schemes for multi- 
dimensional systems of partial differential equations were developed, see e.g. [3,4]. 

The construction of invariant meshless numerical integrators thus seems to be attractive for 
several reasons. From the point of view of invariant numerical schemes, it is beneficial to have 
one more method available that allows one to construct discretization schemes with symmetry 
properties for equations in any space dimension. In turn, from the side of meshless methods 
it is interesting to show that ideas from the field of geometric numerical integration can be 
successfully implemented into such methods. It was shown in this paper that the preservation 
of qualitative properties of a differential equation in a meshless approximation can substantially 
increase the quality of the scheme. The invariant discretization we constructed for a nonlinear 
diffusion equation is able to better reproduce several exact solutions of this equation in practically 
all the parameter ranges that can be tuned in the scheme. We thereby also demonstrated that 
preserving symmetries in a numerical integrator is not solely an academic problem. 

It will be instructive to apply the proposed technique to multi-dimensional discretization 
problems and compare invariant meshless schemes against other types of invariant numerical 
schemes, both in terms of accuracy and computational cost. Also, certain symmetries (e.g. 
Galilean boosts) cannot be preserved on fixed discretization meshes. Moving grid points can 
lead to strongly distorted meshes and are thus on the risk to deteriorate the quality of the 
numerical solution or to slow down the convergence rate. This is what generally happens to 
Lagrangian integration schemes and, as a matter of fact, most invariant numerical schemes pre- 
serving Galilean invariance are Lagrangian integrators [6,8]. On the other hand, it is known that 
certain meshless methods are to some degree insensitive regarding the distribution of the nodes. 
It will therefore be informative to compare invariant meshless methods with discretizations that 
employ classical or invariant moving meshes such as those constructed in [4,11]. 
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